import numpy as np
import matplotlib.pyplot as plt

SE = []

with open("25meV.txt","r") as f:
    for i, line in enumerate(f):
        if i % (192) > 1:
            SE.append([float(x) for x in line.split()])
SE = np.array(SE)
SE = np.reshape(SE,(7,190,3))

def nf(E,T):
    return 1/(np.exp(E/(8.617*10**(-2)*T))-1) + 0.5 +np.sign(E)*0.5

E = SE[0,:,0]
T = np.array([100,200,300,400,500,600,680])


Chi = np.zeros(np.shape(SE))
for i in range(7):
    Chi[i,:,0] = SE[i,:,0]
    Chi[i,:,1] = SE[i,:,1] / nf(E,T[i])



plt.errorbar(SE[0,:,0],SE[0,:,1],yerr=SE[0,:,2],color='navy',marker='o',mfc='none')
plt.errorbar(SE[1,:,0],SE[1,:,1],yerr=SE[1,:,2],color='blue',marker='o',mfc='none')
plt.errorbar(SE[2,:,0],SE[2,:,1],yerr=SE[2,:,2],color='cyan',marker='o',mfc='none')
plt.errorbar(SE[3,:,0],SE[3,:,1],yerr=SE[3,:,2],color='lime',marker='o',mfc='none')
plt.errorbar(SE[4,:,0],SE[4,:,1],yerr=SE[4,:,2],color='gold',marker='o',mfc='none')
plt.errorbar(SE[5,:,0],SE[5,:,1],yerr=SE[5,:,2],color='red',marker='o',mfc='none')
plt.errorbar(SE[6,:,0],SE[6,:,1],yerr=SE[6,:,2],color='darkred',marker='o',mfc='none')


for i in range(7):
    np.savetxt('SE_25meV_'+str(T[i])+'K.txt',SE[i,:,:])

plt.xlim(0,25)
plt.ylim(0,0.002)
plt.legend(["100K","200K","300K","400K","500K","600K","680K"],fontsize=20,frameon=False,ncol=2,handlelength=1.5)
plt.xlabel('Energy (meV)',fontsize=20)
plt.ylabel(r'$S(E)$',fontsize=20)
plt.tick_params(which='both',direction='in',labelsize=20)
plt.savefig('SE.eps',format='eps',bbox_inches='tight')
plt.show()

